This PR consists of several changes that I did in steps.

Step 1

The current approach is integrating the transport velocity through stages. So, between time steps, \tilde{v} is set to v, and then it is integrated for one time step with d\tilde{v}/dt set to the term inside the parentheses here:

Note that this approach is currently flawed, as the dv/dt term in the RHS is missing, so only the physical velocity from the last time step is used at every stage.

In the paper, with the different time integration, this transport velocity is only computed once per step:

While we can interpret this as "transport velocity is integrated within the step starting from v", which leads to the approach explained above, we could also interpret this as "dr/dt = v + δv, where δv is dt * v_c and v_c is the second term inside the parentheses above".
This interpretation leads to a δv = dt * v_c that is either constant throughout the time step and computed between steps, or v_c is recomputed at every stage. The latter is the one I implemented for the following comparisons, but I don't think the difference between computing this only between steps is significant, so we might be able to save some significant computation here.

Results

As can be seen here in the left column, directly behind the cylinder, the old TVF implementation (top row) is less effective at preventing void regions than the new implementation (top row). The old implementation with p_0 = 100k (which is rho_0 * c^2) is about as effective in preventing void regions as the new version with p_0 = 50k, whereas the new version with p_0 = 100k is about as effective as the old version with p_0 = 200k. This indicates that the new implementation might just be a rescaling of the old implementation. We can see here, similar to the results before, that the top left and bottom right show similar particle distribution uniformity (mostly in the red area above the wake, but also in the size of the "reordered" region in front of the cylinder), as well as the top right and bottom left. However, we can now see in the left column that the old TVF implementation changes the solution, producing a slightly wider wake region and lower velocity in the center of the wake. The higher background pressure changes the results even more (bottom right), whereas the top row solution is not affected by the different background pressure. This is a comparison between TVF in this PR (top), consistent shifting from #843 (middle) and inconsistent shifting (bottom).

Step 2

The original 2013 paper was using kernel summation to compute the density. When we use this implementation with the continuity equation, we get the same inconsistency as the inconsistent PST, as particles are now moving with the transport velocity, but the physical velocity is used in the continuity equation.
This can lead to an uncontrolled increase of density in a closed simulation.

In "A generalized transport-velocity formulation for smoothed particle hydrodynamics" (Zhang et al. 2017), the density is integrated with the continuity equation and the transport velocity is used here: